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Abstract 

Hydrodynamic (Landau) instability in combus- 
tion is typically associated with the onset of wrin- 
kling of a flame surface, corresponding to the for- 
mation of steady cellular structures as the stabil- 
ity threshold is crossed. In the context of liquid- 
propellant combustion, such instability has recently 
been shown to occur for critical values of the pres- 
sure sensitivity of the burning rate and the distur- 
bance wavenumber, significantly generalizing previ- 
ous classical results for this problem that assumed 
a constant normal burning rate. Additionally, how- 
ever, a pulsating form of hydrodynamic instability 
has been shown to occur as well, corresponding to 
the onset of temporal oscillations in the location of 
the liquid/gas interface. In the present work, we con- 
sider the realistic influence of a nonzero temperature 
sensitivity in the local burning rate on both types of 
stability thresholds. It is found that for sufficiently 
small values of this parameter, there exists a sta- 
ble range of pressure sensitivities for steady, planar 
burning such that the classical cellular form of hy- 
drodynamic instability and the more recent pulsat- 
ing form of hydrodynamic instability can each occur 
as the corresponding stability threshold is crossed. 
For larger thermal sensitivities, however, the pulsat- 
ing stability boundary evolves into a C-shaped curve 
in the disturbance- wavenumber/ pressure-sensitivity 
plane, indicating loss of stability to pulsating pertur- 
bations for all sufficiently large disturbance wave- 
lengths. It is thus concluded, based on characteris- 
tic parameter values, that an equally likely form of 
hydrodynamic instability in liquid-propellant com- 
bustion is of a nonsteady, long-wave nature, distinct 
from the steady, cellular form originally predicted by 
Landau. 


*This paper is declared a work of the U.S. Government and 
is not subject to copyright protection in the United States. 


1. Introduction 

The notion of hydrodynamic instability in com- 
bustion originated with Landau’s seminal study of 
premixed flame propagation. 1 In that work, it was 
postulated that a flame could be represented by a 
surface of discontinuity, with an appropriate density 
jump across the flame surface, and the normal equa- 
tions of hydrodynamics governing the flow on either 
side of the flame surface. Assuming that this surface 
propagated normal to itself with constant speed, it 
was then determined through a straightforward lin- 
ear stability analysis that a premixed gaseous flame 
was intrinsically unstable to steady (cellular) distur- 
bances whenever, as is generally the case, the density 
on the burned side of the front was less than that of 
the unburned mixture. Since this result is at vari- 
ance with experimental observations of steady, pla- 
nar laboratory flames, it has sparked numerous other 
theoretical investigations that have sought to resolve 
this paradox, through both the inclusion of addi- 
tional effects and the analysis of less phenomenolog- 
ical combustion models. Today, this specific form 
of hydrodynamic instability, generally referred to as 
Landau instability, plays a central role in the study 
of flame/flow interactions. Though Landau’s anal- 
ysis essentially collapsed all thermal effects to the 
change in density across the front, Landau insta- 
bility is nonetheless implicitly associated with the 
thermal expansion of the gas arising from the heat 
release during combustion. Indeed, this instabil- 
ity persists even when the entire flame structure is 
resolved, 2,3 although it is now recognized as a pre- 
dominantly long-wave cellular instability that is sub- 
ject to various stabilizing influences for disturbances 
corresponding to larger wavenumbers. 

A related problem briefly addressed at the end 
of Landau’s original study (for which a more stan- 
dard and less controversial type of result was ob- 
tained), and the one of interest here, was a modest 
modification of the flame-propagation problem just 
described. In this modified problem, the unburned 
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mixture was taken to be a liquid propellant while the 
burned region again consisted of gaseous products. 
The physical existence of a liquid/gas interface led 
to the inclusion of additional physics in the model, 
namely surface tension at the interface and the influ- 
ence of gravitational acceleration (downward propa- 
gation was assumed). As a result, and in contrast to 
strictly gaseous flame propagation, a stability crite- 
rion was derived such that the liquid/gas interface 
was either hydrodynamically stable or unstable in 
the Landau (cellular) sense depending on whether 
the product of the gravitational acceleration and the 
coefficient of surface tension was greater or less than 
a critical value. This result was later extended by 
Levich, 4 who considered the effects of (liquid) viscos- 
ity in lieu of surface tension and obtained a similar 
result for the product of the gravitational acceler- 
ation and the viscosity of the liquid. These classi- 
cal models of liquid-propellant combustion, despite 
their relative simplicity, have long served as seminal 
examples that correctly describe the onset of (cellu- 
lar) hydrodynamic instability in reactive systems. 

Although the assumption of a thin reaction re- 
gion remains a frequently valid and useful simpli- 
fication in many combustion studies, the classical 
assumption of a constant normal burning rate has 
long been regarded as a probable oversimplification 
when applied to the problem of combustion insta- 
bility. Early attempts at modification began with 
the phenomenological assumption of a linear rela- 
tionship between the burning rate and the flame 
curvature, 5 while more modern approaches have em- 
ployed sophisticated asymptotic methods to analyze 
the flame structure and to derive formal expressions 
for the burning rate as a function of local conditions 
at the combustion front, for both gaseous and con- 
densed, combustion waves. 2 ’ 3,6 ' 7 In propellant com- 
bustion, on the other hand, it has long been cus- 
tomary to experimentally measure the pressure re- 
sponse, or pressure sensitivity, of the burning rate, as 
well as (to a lesser extent), its temperature sensitiv- 
ity. Although asymptotic models of propellant com- 
bustion that resolve the combustion-wave structure 
can be developed under various approximations, 
it remains true nonetheless that the direct assump- 
tion of a combustion surface that propagates accord- 
ing to a prescribed burning-rate law, such as a phe- 
nomenological law based on experimental observa- 
tions, allows one to circumvent the intricacies of the 
combustion region and to impose fewer restrictive 
assumptions on the hydrodynamic model in the un- 
burned and burned regions. 

Accordingly, for the case of liquid- propellant 
combustion, the classical Landau/Levich hydrody- 


namic models may be combined and extended to ac- 
count for a dynamic dependence, absent in the orig- 
inal formulations, of the instantaneous burning rate 
on the local pressure and temperature fields. • 
Although exact expressions for the neutral stabil- 
ity boundaries in question can be obtained in spe- 
cific limiting cases, more general results may be ob- 
tained by exploiting the realistic smallness of the 
gas-to-liquid density ratio p. Thus, neglecting ther- 
mal coupling effects, an asymptotic analytical ex- 
pression was derived for the cellular stability bound- 
ary Ap(Jfc), where A p denotes the pressure sensitiv- 
ity of the burning rate and k is the wavenumber 
of the disturbance. 13 The results demonstrate ex- 
plicitly the stabilizing effect of gravity on long-wave 
disturbances, the stabilizing effect of viscosity (both 
liquid and gas) and surface tension on short-wave 
perturbations, and the instability associated with 
intermediate wavenumbers for critical negative val- 
ues of A p . In the limiting case of weak gravity, it 
was shown that hydrodynamic instability in liquid- 
propellant combustion is a long-wave instability phe- 
nomenon, whereas at normal gravity, this instability 
is first manifested through 0(1) wavenumber distur- 
bances. It was also demonstrated that, in general, 
surface tension and the viscosity of both the liquid 
and gas phases each produce comparable stabiliz- 
ing effects in the large-wavenumber regime, thereby 
providing important modifications to the previous 
classical analyses in which one or more of these ef- 
fects was neglected. For A v — 0, the Landau/Levich 
results are recovered from the new model in appro- 
priate limiting cases, although, in the realistic limit 
of small p, this corresponds to a hydrodynamically 
unstable parameter regime. 1 ’ 4,13 

In addition to the classical cellular form of hy- 
drodynamic stability described thus far, there exists 
in this extended model a pulsating form of hydrody- 
namic instability corresponding to the loss of stabil- 
ity of the steady, planar solution to time-dependent 
perturbations. 14 This occurs for negative values of 
the pressure-sensitivity parameter A p , and is there- 
fore absent from the original Landau/Levich models. 
Nonetheless, its existence in the extended model re- 
sults in a stable band of negative pressure sensitivi- 
ties bounded above by the classical Landau type of 
hydrodynamic instability, and below by the pulsat- 
ing form of hydrodynamic instability just described. 
A stable range of negative pressure sensitivities is 
applicable, for example, to certain types of hydroxy- 
lammonium nitrate (HAN)-based liquid propellants 
at low pressures for which nonsteady modes of com- 
bustion have been observed. 15 While it is possible 
that such nonsteady combustion may correspond to 
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secondary and higher-order bifurcations above the 
cellular boundary, 16 it may also be a manifestation 
of the pulsating form of hydrodynamic instability 
just described. 

In the present work, 18,19 we explore the effects 
of incorporating a nonzero temperature sensitivity 
into our previous asymptotic analyses performed in 
the limit of small gas-to-liquid density ratios. 13,14 
This entails a coupling of the energy equation for 
temperature to the previous purely hydrodynamic 
problem, and leads to a significant modification to 
the pulsating boundary such that, for sufficiently 
large values of the temperature-sensitivity parame- 
ter, liquid-propellant combustion can become intrin- 
sically unstable to this newer form of hydrodynamic 
instability. 

2. The Mathematical Model and its Basic Solution 


The extended mathematical model referred to 
above was described previously, 13 but is briefly sum- 
marized here for completeness. In particular, it is 
assumed, as in the classical models, that there is 
no distributed reaction in either the liquid or gas 
phases, but that there exists either a pyrolysis re- 
action or an exothermic decomposition at the liq- 
uid/gas interface that depends on local conditions 
there. For simplicity, it is assumed that within the 
liquid and gas phases separately, the various physi- 
cal properties are constants, with appropriate jumps 
in these quantities across the phase boundary. 

The nondimensional location of the downward- 
propagating liquid/gas interface, which is a function 
of space and time, is denoted by x 3 = $ 5 (xi, x 2 , £)> 
where x 3 is the vertical coordinate and the adopted 
coordinate system is fixed with respect to the sta- 
tionary liquid at i 3 = — oo. Then, in the mov- 
ing coordinate system x = Xi, y = Xi, z ~ x 3 - 
$ s (xx,X2,t), in terms of which the liquid/gas inter- 
face always lies at z = 0, the complete formulation of 
the problem is given by conservation of mass, energy 
and momentum within each phase as 
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where v, © and p denote velocity, temperature and 
pressure, respectively, Pri and Pr g denote the liquid 


and gas-phase Prandtl numbers, p, A and c (used be- 
low) are the gas-to-liquid density, thermal diffusiv- 
ity and heat-capacity ratios, and Fr is the Froude 
number. These nondimensional variables have been 
defined in terms of their dimensional counterparts 
(denoted by tildes) as 
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where U is the reference propagation speed of the 
interface for the case of steady, planar deflagration. 
The nondimensional parameters, some of which first 
appear below in the conditions at the gas/liquid in- 
terface, have been defined as 
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where T u is the unburned (liquid) temperature at 
z = — oo and T a is the adiabatic burned (gas) tem- 
perature at * = -bco, 7 is the coefficient of sur- 
face tension for the liquid surface, and v\ and v g are 
the liquid and gas-phase kinematic viscosities. We 
note that the inverse Froude number Fr~ l repre- 
sents the nondimensional gravitational acceleration 
and that p\Pr g = pPr/, where p = fig/ Pi is the gas- 
to-liquid viscosity ratio. The nondimensional mass 
burning rate appearing in Eq. (9) below is defined by 
Mp\z=o+ i ©iz=o) — j P$ y and 

is assumed to be functionally dependent on the local 
pressure and temperature at the liquid/gas interface. 
By definition, A = 1 for the case of steady, planar 
burning, but perturbations in pressure and/or tem- 
perature result in corresponding perturbations in the 
local mass burning rate. 

Equations (1) - (3) are subject to the boundary 
conditions 


v = 0 , © = 0 at z = — oo , 

0-1 at z = -boo, ©| 2=0 - = ©! 2 =o+ ^ 

and appropriate jump and continuity conditions at 
the liquid/gas interface. The latter consist of con- 
tinuity of the transverse velocity components (no- 
slip), 

n s x v_ = n s x v + , (7) 
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where v± = v| 2=0 ± , conservation of (normal) mass 
flux, 

d<f> m 

n s -(v_ — pv + ) = (1 - p)S($ 3 )-g £- , (8) 

the mass burning rate (pyrolysis) law, 

n s -v_ -S(*.)^f =^(pUo + ,©U=o), (9) 

conservation of flux of the normal and transverse 
components of momentum, 
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and conservation of heat flux 
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where c = c/( 1 - a u ), e is the rate-of-strain tensor 
(e ± = e| 2=0 ±), 7 is the surface-tension coefficient, 
a u is the unburned-to-bumed temperature ratio, and 
5 (<I> S ) and the unit normal n s are defined as 

S(4> 3 ) = [1 + (d$ 3 /dx) 2 + (d^/Sy) 2 ] " 1/2 , (13) 

n s = (-dS./ftr, -d$ a /dy, 1)S($ 3 ) • 
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However, the vector v still denotes the velocity with 
respect to the (x! , x 2 , x 3 ) coordinate system. 

We observe that the thermal and hydrodynamic 
fields are coupled through the temperature depen- 
dence of the mass burning rate A appearing in Eq. 
(9). When A is assumed to depend on pressure only, 
the strictly hydrodynamic problem for p, v and $ 3 
can be analyzed separately [13, 14). In the present 
work, we wish to focus on this more extended cou- 
pling to determine how the hydrodynamic stability 
boundaries are modified when the local burning rate 
depends on temperature as well. In connection with 
this we remark that although the mass burning rate 
A may be typically decomposed as 


4(p|z=o+. Q U=o) 

. rw(i-O(eu =0 -i) 

= A{jp\ z= q+ , ©|*=o) exp ffu + (1 _ ffu) e U= 7j 

(16) 

where A is a rate coefficient and N = E/ R°f a is the 
nondimensional activation energy (P), this more ex- 
plicit representation will not be needed in the stabil- 
ity analysis that follows. Instead, our (linear) stabil- 
ity results may be expressed in terms of the pressure 
and temnerature sensitivities, defined as 


A p — dA/dp\ e=l 

-4© = dAj 9 ©| @ _ lp _ 0 = ~ “b ’ 

respectively. However, we note that if the nondi- 
mensional activation energy is large, then the last of 
Eqs. (17) implies that A e would likely be larger in 
magnitude than a fact that will play a role in the 
relative scalings introduced for A p and A© below. 

A nontrivial basic solution to the above prob- 
lem, corresponding to the special case of a steady, 
planar deflagration, is given by 


Here, the factor multiplying 7 in Eq. (10) is the 
curvature -V ■ n s of the liquid/gas interface in the 
moving coordinate system, and the corresponding 
expressions for both the gradient operator V and 
the Laplacian V 2 in this system are given by 
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The remainder of the paper is devoted to a linear 
stability analysis of this solution. 

3. The Complete Linear Stability Problem 

Before proceeding with any further approxima- 
tions, we determine in a straightforward fashion the 
linear stability problem and the corresponding dis- 
persion relation (in the form of a system of alge- 
braic equations) that govern the behavior of har- 
monic perturbations about the basic solution (18). 
In particular, in terms of the perturbation quanti- 
ties (f> s {x , y,f) = ^ s (t), u(x,y, 2 ,f) ~ 

v(x,y,z,t) <( x,y,z,t ) = p(x,y, z,t) - p°(z) 

and 0(x,y, z,t) — © — ©°(z) — (f) s dQ°/dz , the prob- 
lem obtained when Eqs. (1) - (3), (6) - (12) are 
linearized about the basic solution (18) is given by 
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where c = c(l - u u ) 1 . 

Nontrivial harmonic solutions for <j> s , u and 
proportional to v, that satisfy Eqs. (19) 

- (21) and the boundary/boundedness conditions at 
z = ioo are given by 
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for z > 0, where we have normalized the above so- 
lution by setting the coefficient of the harmonic de- 
pendence of <t> s to unity. Here, the signs of k\ and 
k 2 may be either positive or negative and we have 
employed the definition k = y/kj T k 2 , and p, q } r 
and s are defined as 

2p = 1 + y/l + 4(*w + k 2 ) , (37) 
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2Pri q= 1 + \A + 4Pn(ia; + Pri fc 2 ) , (38) 

2 uPri r = 1 - \/T+ 4/LiPrj (iw/j + /iPri A: 2 ) , (39) 

2pAs = 1 - a/1 4- 4p 2 A(i<u + Afc 2 ) , (40) 

where we have used the fact, noted below Eq. (5), 
that pAPr 9 = /iPrj . 

Substituting this solution into the interface con- 
ditions (23) - (30) and using Eq. (19) for z > 0 yields 
eleven conditions for the ten coefficients i> i - b io and 
the complex frequency (dispersion relation) iu>(k). 
In particular, these conditions are given by 
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While the above problem is linear in the coefficients 
61 - 6a, the relationship for iu is highly nonlinear. 
Accordingly, we seek asymptotic solutions for the 
neutral stability boundaries in the limit that the 
gas-to-liquid density and viscosity ratios p and p are 
small, as is Fr~ l in the case of reduced gravity. We 
shall also introduce appropriate scalings for the sen- 
sitivity parameters A p and A©, where the limiting 
results for A© = 0 were derived previously [13, 14]. 
In particular, we introduce a bookkeeping parameter 
e«l and consider the realistic parameter regime 


p — P* e , M = Fr 1 


g ~ 0(1) 
g*e ~ 0(e) , 


(52) 


where the the upper (lower) choice for Fr~ l corre- 
sponds to normal (reduced) gravity We note that 
based on characteristic densities for liquid propel- 
lants and product gases at elevated temperatures, a 
realistic value for e is on the order of 10” 3 or 10 4 . 
In this parameter regime, it was seen that the ap- 
propriate scaling for A p to describe the neutral sta- 
bility region was A p = A*e, whereas the appropri- 
ate scale that describes the fully-developed effects of 
thermal coupling turns out to be A© = A@€ lyM , as 
shown in Section 6 below. Thus, we observe that 
A e /A p r 0(e -3 / 4 ) > 175, which is at least several 
times as large as typical values of the nondimen- 
sional activation energy N, but still consistent with 
the last of Eqs. (17). However, as shown in Section 
7, the first significant effects of thermal coupling are 
in fact felt on an intermediate scale Aq ~ 0(e x / 2 ), 
in which case Aq/A v ~ 0(e -1 / 2 ) > 30, which is 
a quite typical value for N_. Thus, the parame- 
ter regime of interest for describing modifications to 
the hydrodynamic stability boundaries arising from 
temperature-sensitivity effects coincides with realis- 
tic values of the sensitivity coefficients. 


4. The Inviscid Limit 

Although both liquid and gas-phase viscous ef- 
fects were shown to be generally comparable to those 
due to surface tension when thermal coupling was 
neglected, 13,14 the qualitative nature of the cellu- 
lar boundary remains preserved in the zero- viscosity 
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limit, as does the existence of the pulsating bound- 
ary. Accordingly, we shall, for simplicity, mainly 
consider the effects of thermal sensitivity on hydro- 
dynamic stability in the inviscid limit, reserving the 
detailed consideration of thermal coupling in the al- 
gebraically more involved viscous case for a future 
study. 19 In particular, before introducing the scal- 
ings indicated above, we first set Pn = PTg — 0* 
Consequently, since q — > oo according to Eq. (38), 
we set 63 = 65 = 67 = 0 in Eqs. (33) - (35), and 
observing that Eqs. (48) - (49) then reduce to the 
no-slip conditions (43) - (44), we thus eliminate Eqs. 
(41), (48) and (49) from the above set. Using the 
fact that r — * — iup , the solution for the remaining 
coefficients in terms of iu is given by 

b * 

8 2 p(iup - k)(iup 4- fc) 

-{(i -p) 2 [MV + * 2 ] (53) 

+ p( 1 - p 2 )Fr~ 1 k + p(l + p)7fe 3 |, 


61 = - ,. W + ^ 2 - [(iu) + k)(l - p)k - p{iup - k)b$] , 

(54) 

h = - if * K*"'’ - t)a - rtt + A <* + *)>»] • 

H ' +f,)k (55) 

(56) 


, vpki L L _ upk 2 u 
O 4 = yr-O g , ^6 = , 


A©6l0 = Ar 


it 2 
pFr~ l 


k 2 

- fc) 2 (l - p) 


+ 


&9 = &10 “ 1 + 


P( 1 + 

p(ia?p - fc)(iu> 4- k) 

+ (14 p)fc 2 

(1 - p)k — 2ia;p p(io;p - fc) 


^8 


14 -p 


(1 4- p)k 


h, 


(iu 4" fc)(l — p)k — p(iup — fc)6g 


(57) 


(1 4- p)(icj - fc)fc 


(58) 


Substitution of these results into Eq. (51) then gives 
a single implicit equation for iu given by 


[c(pAs — 1) 4- 1 - p]b\ 
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iu 4- k 


4- 


1 -bp 
P 

(1 + p)k 


iu — k 
iup - k 


iu — k 


(jp-k- 1) - (1 + c + cp)iu - k 
u 

(p - k - 1) - ( iup - k) 


= 1 —p + iu [l + 6(1 - p)] , 

(59) 

where 6 8 (zw) and bi 0 {iu) are given explicitly by Eqs. 
(53) and (57), provided Aq ± 0. Equation (59) is 


thus the dispersion relation that determines the neu- 
tral stability boundaries in that case. In the limit 
Aq -* 0, the dispersion relation is given by setting 
the right-hand side of Eq. (57) to zero, in which 
case Eqs. (58) and (59), which are then decoupled 
from Eqs. (53) - (57), determine 69 and &io- Al- 
ternatively, we may solve Eq. (59) for b\o(iu)*&nd 
substitute this result into Eq. (57), which then be- 
comes the dispersion relation for arbitrary A©. We 
remark that the special case A© = 0 yields an an- 
alytically tractable solution for the neutral stabil- 
ity boundaries, 12 but since all other results to be 
presented are obtained for the parameter regime de- 
fined by Eqs. (52), we shall restrict ourselves to this 
regime for this limiting case as well. 

5. Hydrodynamic Stability Boundaries 
in the Limit A© = 0 


In the limit A© = 0, the results obtained from 
the dispersion relation (57) in the parameter regime 
(52) are as follows. 13 First (see Figure 1), the pul- 
sating hydrodynamic stability boundary is given, to 
leading order in e, by A* = — p* for all disturbance 
wavenumbers [the exact relation is given by A* = 
-p*/(l - ep*)], with instability occurring below this 
critical value. For the cellular stability boundary, 
the corresponding leading-order expression depends 
on the magnitude of the wavenumber k . In partic- 
ular, there are three wavenumber scales to be con- 
sidered; the 0(1), or outer, scale fc, a far outer scale 
kf = fee, and an inner scale defined by either fc* = 
k/e if Fr- 1 ~ 0(1) or k t = k/t 2 if Fr~ l ~ O(e). 
In each of these regions, the cellular stability bound- 
ary, on which the complex frequency iu is identically 
zero, is given by 



A;~A; U) ~± P *(p*ik f -i), (6i) 

and 



/ P*(P*9 ~ ki)/2ku 
\ P*(P*9* — ki)/2ki, 


Fr' 1 ~ 0(1) 
Fr- 1 ~ 0(e) , 
(62) 


respectively. Matching the above solutions to one 
another, a uniformly valid composite expansion, de- 
noted by A*( c \k), can be constructed in the region 
A* < 0 (the basic solution can be shown to be un- 
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stable for A* p > 0). In particular, 



a; {1) + a; {0) + a; U) 

- lim - lim A * (/) 

fci-. oo P fc/-*0 p 


~ ~\p* + £«P*V 

, / ep* 2 ff/2A:, Fr' 1 ~ 0(1) 

+ \ e 2 p* 2 g* /2k, Fr 1 ~ 0(e) , 


(63) 


where the definitions of ki and kf have been used 
to express the final result in terms of k. We note 
that this asymptotic result, exhibited in Figure 1, 
is far simpler than the exact result, 12 which, in our 
present notation, is given in the region A* < 0 by 


a;/ p * - 

p»c( 1 - p*e)Fr~ l + p*ejk 2 - (1 - p*e)k 

p* 2 e 2 (3-p*e)Fr- l +p* 2 e 2 ik 2 + (l-p*«X2-p*c)fc 

(64) 

It is thus easily seen from Figure 1 that there 
is a stable region between the pulsating and cellu- 
lar stability boundaries for A© = 0, and this re- 
sult is preserved (in fact, enhanced) when viscous 
effects, which play a stabilizing role with respect 
to both boundaries for sufficiently large wavenum- 
bers, are included in the analysis. 13 - 14 We observe 
that for the Landau (cellular) form of hydrodynamic 
instability, gravity and surface tension are stabiliz- 
ing with respect to sufficiently small and sufficiently 
large wavenumber disturbances, respectively. The 
essential difference, as illustrated in Figure 1, be- 
tween the normal and reduced gravity limits in the 
realistic parameter regime considered here is that in 
the latter instance, gravity is only capable of sta- 
bilizing disturbances whose wavenumbers are very 
small, 0(e 2 ), whereas in the former case, gravity 
is sufficiently strong to stabilize disturbances whose 
wavenumbers are 0(e). As a consequence, hydro- 
dynamic instability becomes a long-wave instability 
phenomenon in the reduced gravity regime consid- 
ered here, since, from Eq. (63), the most unstable 
wavenumbers are 0(>/c), rather than 0(1), in that 
case. We remark that although positive values of the 
pressure sensitivity A p are hydrodynamically unsta- 
ble, zero and negative values of A p over certain pres- 
sure ranges are, in fact, characteristic of the so-called 
“plateau” and “mesa” types of solid propellants, 17 
as well as for the HAN-based liquid propellants men- 
tioned above. 15 



Figure 1. Inviscid hydrodynamic stability boundaries 
in the limit of zero thermal coupling (Aq = 0 / 

6. Hydrodynamic Stability Boundaries for Aq > 0 

Retaining the parameter scalings introduced in 
Section 3, we first consider the 0(1) wavenumber 
regime. For the case A© = 0 considered previously, 
it turns out 14 that the appropriate expansion of the 
dispersion relation iu(k) is of the form 

iuj ~ t~ l f 2 (i<jjQ *f e 1 ^ 4 ^! e l f 2 iu)2 + ■*■), (fib) 

where the leading-order term was suggested by the 
explicit results available for the inviscid case, and 
the expansion in quarter powers of e is suggested 
by the leading-order term in the resulting expansion 
for p given below. We now introduce the scaling 
Aq = e n A* e , where the value of n > 0 is to be deter- 
mined so as to give a nontrivial modification in the 
dispersion relation at the appropriate order, corre- 
sponding to the first effects of a nonzero temperature 
sensitivity on the neutral stability boundary for 0(1) 
wavenumbers. 
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Substituting Eq. (65) and the various scalings 
into Eqs. (37), (40) and (53), expanded forms for p, 
s and b$ are obtained as 

p ~ Poe _1/4 +Pi+P2« 1/4 + ■ ■ , s~0(e 1/2 ), (66) 

(iw 0 ) 2 p* + k 2 _i 1 _ 3/4 

2p*k € -^(wv 0 )^i)e ' 

~ ^ [( iw i) + 2{iu>o)(iuj 2 )] e~ 1/2 + • • • , 

(67) 

where 

Po = (*wo) 1/2 , Pi = | [fwi/(iw 0 ) 1/2 + lj , 

P2 = ^(^o) -1/2 [l + 4 k 2 + 4iw 2 - (iwi) 2 /iw 0 ] • 

(68) 

Consequently, from Eqs. (57) and (67), 6 10 has the 
expanded form 

b io ~ 

+ ~ [(two)V - fc 3 ] - j-(iuo)(iw,)e 1/4 

+ §[^{MV-3* 2 } 

— ( iuf \ ) 2 — 2(icJo)(^2)j + * * * 

(69) 

which, when substituted into the dispersion relation 
(59) along with the other expansions deduced thus 
far, yields a sequence of equations for the deter- 
mination of the coefficients in the expansion (65). 
Thus, at leading order, 0(e~ n ~' 1 ^ 4 ) 1 Eq. (59) yields 

/n\ 


b\o = 0, which, from the implied definition of b\V 

in Eq. (69) above, gives 



(70a) 

° r 2 


(iojo) 2 = (2^4p + P # ) • 

(706) 


Hence, iu)Q is identically zero for A* = — p*/ 2, which 
corresponds to the leading-order hydrodynamic cel- 
lular boundary (60) in the 0(1) wavenumber regime 
when Aq = 0. Values of A? > —p*/ 2 imply cel- 
lular instability, since in that case there exists, ac- 
cording to Eq. (70b), a root of the dispersion re- 
lation such that TZe{iuJo} > 0, Xm{i<jJo} = 0. On 


the other hand, the stability of the region A* < 
—p* /2 is still indeterminate, since Eq. (70b) in- 
dicates that 7£e{zwo} = 0 there. Hence, the next 
nontrivial term in the expansion of iu> is required to 
determine whether this region is stable or unstable, 
although the fact that Xm {zuo} ^ 0 implies that 
disturbances have a pulsating character for values of 
A p below the cellular stability boundary. 

At the next order, 0(e _n ), in the analysis of 
the dispersion relation, the equation 6$ = 0 is ob- 
tained, whence the definition of b[ X Q in Eq. (69) im- 
plies that iu i = 0. Proceeding with the next higher 
order equation obtained at 0(e“ n+1/4 ) from the ex- 
panded version of Eq. (59), it can be seen that terms 
other than those proportional to A^~ l will appear 
provided the choice n = 1/4 is made. In that case, 
collecting terms of 0(e°) and using the previous re- 
sults = 0 leads to the equation 

-Po&m — 2fc ^ Z ] = 0 j (71a) 

which, upon using the definitions of po and 6^, may 
be solved for zaJ 2 as 

i(jJ2 = 

[(«■*> ) V - 3 k 2 ] [A; + p* + fci4@(iu;o) -3/2 ] 

• [a; + P * + P * 3/2 k- 1/2 A' e (2 a; + p*) _3/4 ] , 

(716) 

where the expression (70b) for iu 0 has been used 
to obtain the final equality. It is readily seen that 
for Aq — 0, 7^e{zo;2} - 0 for A* = -p* < -p*j 2. 
.That is, as described in the previous section, there 
exists, since Xm{zu;} ^ 0, a pulsating neutral sta- 
bility boundary that lies below the cellular bound- 
ary such that the region between the two, namely 
-p*/2 < A* < -p*, is stable [7Ze{iu2] < 0], and 
the region below the pulsating boundary is unstable 
[7£e{up2} > 0] . 

For Aq > 0, the pulsating boundary is modi- 
fied. In particular, in the region A* < -p* /2 be- 
low the cellular boundary, the principal value of the 
complex factor in Eq. (71b) may be expressed as 
(2A* p + p*)-V 4 = [~(2A* + p*)] -3/4 e -3 ’*/ 4 , and 
thus the neutral stability condition 'Jlefyuz} — 0 is 
given by 

A;+p’-^V2p' 3/2 k- l / 2 A’ e [-( 2 a; + p *)] - 3/4 = 0 . 

( 72 ) 
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Writing A * = -p* 4- ap* / 2, where a represents the 
deviation, in units of p*/ 2, in the pulsating sta- 
bility boundary from its value in the zero-thermal- 
sensitivity limit Aq = 0, the condition (72) becomes 
a - y/2p*~ l / 4 k~ l / 2 AQ(l - a)~ 3 / 4 = 0, or equiva- 
lently, 

a 4 (l -a) 3 = V _1 fc- z >l^ 4 . (73a) 

Equation (73a) is an implicit expression for the pul- 
sating stability boundary a(fc), where the value a = 
1 corresponds to the cellular boundary and a = 0 is 
the pulsating boundary for A$ = 0. Alternatively, 
writing A* = {—p*j 2)(1 + 6), where 6 = 1 — d is the 
negative deviation, in units of p* j 2, from the cellu- 
lar stability boundary A* = — p*/2, Eq. (73a) can 
be rewritten in terms of b as 

S 3 (l -6) 4 =4p* _1 fc- 2 ^ 4 , (736) 

where b = 1 (6 = 0) corresponds to the pulsating 
(cellular) boundary in the limit Aq = 0. A plot 
of fc(6) for several different values of A* e is shown in 
Figure 2, which, when rotated -90° so that the k - axis 
is horizontal, is readily interpreted in the context 
ofPigure 1, the lines A* = — p*/ 2 and A* = — p* 
in that figure corresponding to b = 0 and 6 = 1, 
respectively, in Figure 2. 

Referring to Figure 2, it is seen that for A@ > 
0, the pulsating boundary becomes C-shaped (in 
the rotated frame of reference), the upper branch 
approaching the cellular boundary b = 0 as k — > 
oo, and the lower branch approaching the original 
(A@ = 0) pulsating boundary 6=1. These same 
limits are approached for any fixed value of A: as 
A@ — ► 0. The portion within the C-shaped curve 
is the stable region, and thus ..not only is steady, 
planar burning intrinsically unstable for sufficiently 
small wavenumbers, but, for finite values of k } any 
crossing of the C-shaped boundary from the stable 
to the unstable region corresponds to the onset of 
a pulsating instability. As A@ increases, the turn- 
ing point of the C-shaped pulsating boundary, which 
occurs at b = b c = 3/7 (i.e., at A* = -5p*/7) cor- 
responding to the critical wavenumber k = k c = 
(343/72)(21/p*) 1 '' 2 J 4£ 2 « 21.83 ^q 2 /p* 1/2 , shifts to 
larger values of k as A@ increases. On the other 
hand, as A@ becomes small, the turning point shifts 
to small values of k such that k c eventually leaves 
the 0(1) wavenumber region for which Eqs. (73) 
are valid. Thus, as A@ becomes small, the original 
pulsating and cellular boundaries are recovered in 
the 0(1) wavenumber regime, but as Aq becomes 
large, the original cellular boundary lies within the 


unstable region for 0(1) wavenumbers, and the basic 
solution becomes intrinsically unstable to oscillatory 
disturbances.* 


Pulsating Stability Boundary 
Ae - 0(e IM ). k - 0(1) 



Figure 2. Pulsating hydrodynamic stability bound- 
aries for k ~ 0(1) and Ae ~ 0(e lj/4 ). 

A composite asymptotic solution for the neu- 
tral stability boundary in the regime Aq ~ O^ 1 / 4 ) 
is thus obtained by matching the cellular and pulsat- 
ing boundaries in the far outer wavenumber regime, 
where the former is given by Eq. (61) and the latter 
is given trivially by A* = -p*, with the appropriate 
solution branch of Eq. (72) in the 0(1) wavenum- 
ber region. In particular, we denote the two so- 

*These results correct an erroneous conclusion drawn 
in a previous study, 12 where an algebraic error in a 
similar, but less formal, calculation incorrectly sug- 
gested that the pulsating boundary recedes (to more 
negative values of A*) as A 0 increases, leaving the 
cellular boundary unaffected. 
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lution branches of Eq. (72), which correspond to 
the portions of Figure 2 that lie to the left and to 
the right of the minimum at b = b Ct by A*^° ,u) (A;) 

and A*(°' l \k), where the superscript V’ denotes, 
as before, the outer, or 0(1), wavenumber region 
and the superscripts “u” and “1” denote the up- 
per and lower (rotate Figure 2 by -90°) portions of 
the double-valued pulsating boundary A*(k). These 
branches are given implicitly by Eq. (73b), where, 
by definition, 6 = — (2A*/p* -f 1). Along the upper 

branch, A*^ 0,t ^ — > —p*f 2 (z.e., b — > 0) as k — oo, 
which can be matched with Eq. (61) since, from 
Eq. (61), A*^ — * —p*/ 2 as kf — ► 0. Similarly, 

_p * (j . e .,6 —> 1) as k — » oo, which clearly 
matches the pulsating boundary A* = —p* in the 
far outer wavenumber region. As a result, a leading- 
order composite stability boundary spanning both 
the outer and far outer wavenumber regions is given 
by 


Kik) 


/i;^>(fc) + e p* 2 7 fc/2. 


a; > ~s P */7 

A;<-5p*/7, 

(74) 


Figure 3. Composite hydrodynamic stability bound- 
ary for A© 0{e^ 4 ). 


where the second term in the top expression has been 
expressed in terms of the outer wavenumber variable 
k. This composite boundary is shown in Figure 3. 
We note, based on the above construction, that the 
lower branch of Eq. (74) is a pulsating boundary for 
all wavenumbers, whereas the upper branch transi- 
tions from a pulsating boundary for 0(1) wavenum- 
bers to a cellular boundary for 0(€ _1 ) wavenumbers. 
Indeed, from Eq. (70b), we observe that in the outer 
wavenumber regime, the size of the upper region of 
oscillatory instability, which is bounded below by 
the upper branch of the pulsating stability boundary 
and above by the region of nonoscillatory instability 
beyond the old cellular boundary A * ~ -p*/ 2 for 
Aq = 0, shrinks to zero as k becomes large. 

7. Evolution of the Pulsating B oundary for Aq 1 

The nature of the evolution, as Aq decreases, of 
the pulsating stability boundary depicted in Figure 
2 to that shown in Figure 1 for the case Aq == 0 
may be determined by analyzing the dispersion re- 
lation for smaller order-of-magnitude wavenumbers 
and appropriately rescaled values of A©. In par- 
ticular, the results (70b) and (73b) suggest that if 
we consider k ~ 0(e°), then A@ ~ 0(e a/2 ) and 
zc^o ~ 0(e a ), where the latter imply that A© and iu> 


Composite Pulsating/ Cellular Stability Boundary 
Aq = .5, p* = y = l t e - .005 
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are 0(e< 2a+1 ^ 4 ) and 0(e (2a “ l) / 2 ), respectively. A 
preliminary analysis then suggests that qualitative 
deviations from the C-shaped nature of the pulsating 
boundary exhibited in Figure 2 occur for wavenum- 
bers of order e 1/2 , which is another intermediate 
wavenumber scale that lies between the inner and 
outer wavenumber regions that were defined in. the 
previous section. The above analysis is thus re- 
peated for k — ke 1 / 2 and A© = A©e 1/f2 , the other 
parameters remaining unchanged from their previ- 
ously assumed orders of magnitude. Consequently, 
in place of Eq. (65), a solution for the dispersion 
relation is now sought in the form 

ioj ~ id ) o 4- e l ^ 4 id ) i 4 e l ^ 2 id ) 2 4 . (75) 

Substituting Eq. (75) and the new scalings in- 
troduced above into Eqs. (37), (40) and (53), the 
new expanded forms of p, s and are given by 

p ~ po 4- pie 1/4 4- p 2 e 1/2 H 1 s ~ 0(c) , (76) 

(id)o) 2 p* + k 2 _i /2 1 w.~ \ - 1/4 

2p*k k 

f(tcii) 4 2 (iwo) (i^2 )] €° 4 * • • , 

2k , , 

(77) 

where 

p° = ~ [l 4 (zu>o) 1/2 ] 1 Pi — + 4iLj 0 )~ 1 / 2 , 

p 2 = + 4^o) _1/2 - 

(78) 

As a result, Eqs. (57) and (67) give the correspond- 
ing expansion of 610 as 


610 ~ 

e°S (4 + ** [(**,) V - 3* 2 1 ^ + • • • 

+ ~ TT f(-o)V - k 2 } - 4-(tw 0 )(iwi)e l/4 

A© l 2k L J k 

+ 4[^{ (i( i 0 )V-3ik 2 }-(«D 1 ) 2 
2 fc ^ k 

- 2(t(2»o)(jw2)je 1/2 H | 

~^ o 4S1 1 0 V/ 4 + 6< 2 o V/2 + ... , 

(79) 

Substituting these expansions into the dispersion re- 
lation (59), one obtains from the leading 0(c°) terms 
in that equation the result = 0, or, according to 
Eq. (79), 

(^o) 2 = ^ (2A; + P*) . (80) 


Comparing Eq. (80) with Eq. (70b), it is thus seen 
that the same leading-order result is obtained for 
iu 0 (k) for Ae ~ 0(e 1/r2 ) as was obtained for i uo{k) 
for Aq ~ 0(e 1 / 4 ). In particular, id)o is real and 
positive for A* > —p*/ 2, and is pure imaginary for 
A* < — p*/2. Thus, the prediction of cellular insta- 
bility is recovered in the former case, but stability 
in the latter region can only be ascertained by cal- 
culating the next nontrivial term in the expansion 
(75). 

Collecting terms of 0(e 1/?4 ) in the expanded ver- 
sion of Eq, (59), it is concluded that b$ = 0, which, 
from Eq. (79), implies that %Cj\ = 0. Proceeding to 
the next order, we collect terms of 0(e 1//r2 ) in Eq. 
(59), which gives rise to the nontrivial equation 

(1 - c - po) + (2k + p*6f) (l + = 0 , 

(81a) 

where, from Eq. (77), 6g 0 ^ = — {idj^){id) 2 )/k. Upon 
substituting the expressions for and (iuo) 2 from 
Eqs. (79) and (80), we obtain the expression for id ) 2 
as 


id)2 = k(Ap — p*) j^A* 4 p* 

Aep * 2 { 14 2id)o — \/l 4 4zc2?o ^ 

k(2A* 4 p*) V 1 - 2c - y/1 4 iid) 0 ) 

(816) 

where, from Eq. (80), za>o = (k/p*)(2A* 4 p*) 1/2 = 

{kfp*)[ - (2 A* 4 p*)] l/2 e i * /2 ) the second equality 
denoting the principal root when A* < -p* / 2. The 
expressions (81b) and (71b) for the two growth-rate 
corrections id )2 and iu 2 in their respective wavenum- 
* ber regimes collapse to the same result in the limit of 
zero thermal sensitivity, yielding the same pulsating 
stability boundary A * = — p* as before. However, 
for nonzero A©, the forms of these two results differ, 
leading to qualitative differences in the correspond- 
ing neutral stability boundaries. In particular, it will 
be shown below that the pulsating neutral stability 
boundary obtained from Eq. (81b) approaches that 
shown in Figure 1 (i.e., A* = - p *) in the limit of 
zero thermal sensitivity, whereas it approaches the 
form shown in Figure 2 as A© becomes large. That 
is, the evolution between these two limiting forms 
occurs on the 0(e 1 ^ 2 ) wavenumber scale for 0(c x / 2 ) 
values of the thermal sensitivity parameter A© . 

Setting the real part of iu 2 to zero, the equation 
governing the (pulsating) neutral stability boundary 
A*(fc) in the present wavenumber regime is deter- 
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mined from Eq. (81b) as 


by 


K+p' 

A&p* 2 ( 1 4- 2iwo — y/1 -f- 4 z&q 1 _ q 

k{2Ap 4- p*) 1 1 — 2 c — vTT 4ia> 0 / 

( 82 ) 

which, since i(b o depends on A * and is a highly 
nonlinear equation for A* v (k). Thus, in order to de- 
termine the nature of the stability boundary as a 
function of A©, it is useful to consider certain lim- 
iting cases. We first consider possible solutions of 
Eq. (82) in the large (scaled) wavenumber regime 
k » 1. In that limit, we observe from Eq. (70b) 
that iu>Q ~ O (k) » 1, and thus a first approxima- 
tion to Eq. (81) in this regime is given by 


K + p' 


AeP * 2 
k(2A; + p*) 




(83a) 


or, using the expression for the principal root ta>o 
given below Eq. (81b), 

A;+p*-±V2p' 3 ' 2 k-'/*Ae [-(2 a; + p*)]~ 3/4 ~ 0 , 

(836) 

which, since k~ l / 2 A& = {k / t l f 2 )~ l f 2 [A® / t l J 2 ) = 
k~ l/2 {A e /e x / 4 ) = k- l ? 2 A%, is identical to Eq. (72). 
According to the discussion of Eq. (72) above, the 
turning point in the corresponding neutral stability 
boundary exhibited in Figure 2 occurs at the value 
k c oc A 0 2 . Thus, we conclude that Eq. (83b) is a 
valid approximation for Eq. (82) in the large (scaled) 
wavenumber regime k » 1 provided A| > 1, where 
the latter is a necessary condition to insure that k 
is large everywhere along the stability boundary de- 
scribed by Eq. (83b). For such values of A©, the 
pulsating stability boundary for large k approaches 
that illustrated in Figure 2 (with k replaced by k 
and A 0 replaced with A©). 

In addition to the boundary described by Eq. 
(83b) and Figure 2 for large k and A©, there is 
another pulsating branch that is present for small 
wavenumbers, as suggested by the observation that 
there exists a solution of Eq. (82) such that A* — > 

—p* as k — + 0, irrespective of the magnitude of A©. 
To simplify the analysis and further discussion, we 
mainly consider a simplified form of Eq. (82) by 
restricting further consideration to the reasonable 
limit in which the heat-capacity ratio c is small. 
Since 1 4- 2io)o — \/l 4- 4fa)o = (1 — \/l""-f- 4icI>o ) 2 /2, a 
first approximation to Eq. (82) in that limit is given 


A; + p ' + i W a;+ p -) ^ + 4 “° ' 1 ) ~ 0 ' 

c -+ 0, 

t . (84) 

which clearly gives the same limiting approximation 
(83a) for large k . To analyze Eq. (84) for arbitrary 
k t it is convenient to convert it to an alternative 
form as follows. Writing 1 + 4ztI>o in the polar form 
1 4- 4iw 0 - r cost? 4- isint?, where r 2 = 1 4 - 16cjq, 
cost? = 1/r and sint? = 4cj 0 /r, Eq. (84) may be 
rewritten as 


2u 2 0 (A; + p*) = kA Q 



(i + 1 6w%y/ 2 + 1 
2 



(85) 

where, from Eq. (80), = -(k/p*) 2 (2A* + p*). 

Substituting this expression for w 2 into Eq. (85) 
and rearranging/squaring the latter twice so as to 


remove the fractional powers, a polynomial equation 
for the inverse relation k(A p \ A e ) is obtained as 


4 k 2 (2A* p + p*) 3 (A* p + p*) 4 

-s fcVMe {2A; + p*) 2 (A; + p') 3 
+ k [5 p* 4 i| (2 a; + p*) (a; + p*f + p* 6 i|] 
-p* 6 i|(x; + p*) =o, 

(86a) 

or, introducing b = ~{2A* p j p* + 1) as before, 


P6 3 (l - 6) 4 -Mfc 2 i e 6 2 (l - 6) 3 

+ k [5i|6(l - b) 2 - 4X|/p*] + 2X|(1 - 6) = 0. 

1 (86 b) 

An even more compact version is obtained by defin- 
ing k = 6(1 - 6)fc/A©, giving the cubic equation 

P 4- 4p + (5 - a)k 4- 2 = 0 , (86c) 

where a(6) = 4A 0 /[p*6(l - 6) 2 ]. However, since it 
is the relation A*(&;A©), or equivalently, k(b; A©), 
that is ultimately desired, we shall mainly consider 
the form given by Eq. (86b). 

In the regime A© 1 there exists a solution 
k ~ O(A 0 ) as suggested above. More precisely, this 
solution of Eq. (86b) may be sought in the form 

k ~ A|(fc 0 + Aq 1 k\ 4- A e 2 k 2 4" • ) . (87) 


whence, upon substitution of this expansion into Eq. 
(86b), the coefficients k t are determined recursively 
as 


h = 26- 3 / 2 (l -6)- 2 /p* 1/2 , ki = — 26 -1 (l - b )- 1 , 

(88) 
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and so forth. As expected, the leading-order approx- 
imation k ~ AqIcq is identical in form to the result 
(73b) and the next term in the above expansion thus 
provides a negative correction to this result. On the 
other hand, for large A© there also exists a solution 
k 1 of the form 

k ~ "b A^k\ + H ) ■ (89) 

Substituting the latter into Eq. (86b), the ki are 
also determined recursively, the first few of which 
are given as 

k 0 =p*(l-b)/2, k\ — 0 , k 2 = 5p* 2 S(l - 6) 3 /8 , 

(go) 

where the last of these provides, for 0 < b < 1, a 
strictly positive correction to the leading-order ap- 
proximation k (2 A©) l p*( 1 - 6). Although the 

branches described by Eqs. (87) and (89) were ob- 
tained for A© > 1, they remain reasonably accurate 
representations of the actual branches even for 0(1) 
values of this parameter, as shown in Figures 4a, b 
(discussed immediately below). 

As A© decreases through 0(1) values, the turn- 
ing point in the stability boundary for large k occurs 
for smaller values of fc, and eventually the expansion 
for large k ceases to be a valid approximation of Eq. 
(86b). Simultaneously, the branch corresponding to 
small values of k reaches larger values of k such that 
the above expansion for small k ceases to be valid 
as well. These two branches eventually intersect, 
restoring a region of stability for a range of negative 
values of A*. This development is illustrated in Fig- 
ures 4a-f, which were obtained directly from the ex- 
act roots of the cubic equation (86b) with respect to 
k [for comparison, we also indicate the asymptotic 
representation for the upper branch corresponding 
to one- and two-term truncations of Eq. (87)]. We 
observe that after the two pulsating branches inter- 
sect (as, for example, in Figures 4c-e), there exist 
two critical values of k corresponding to the onset of 
pulsating instability at corresponding critical values 
of b y or equivalently, A*. That is, in this range of 
temperature sensitivities, a pulsating hydrodynamic 
instability occurs for either sufficiently small or suf- 
ficiently large negative values of A* as either the left 
or right (upper or lower in the rotated frame of refer- 
ence) stability boundary is crossed. As Aq decreases 
further, the separated branches move further apart 
until ultimately, in the limit Aq — + 0, we recover 
from Eqs. (86) the two roots 6 = 0 and 6 = 1 corre- 
sponding to the cellular (A* = — p*/ 2) and pulsat- 
ing (A* = — p*) stability boundaries in the present 


[k ~ 0(e 1/2 )] wavenumber regime. Thus, as the 
temperature-sensitivity parameter A© decreases to 
small values on the scale A©/A p 0(e~V 2 ) 50, 

the region of hydrodynamic stability depicted in Fig- 
ure 1 is recovered. 


Pulsating Stability Boundary 



Figures 4a-f. Pulsating hydrodynamic stability 
boundaries for k = fce 1 / 2 ^ 0(e 1 ^ 2 ) and decreas- 
ing values of A© = A©e l ^ 2 ~ 0(e 1 ^ 2 ). The broken 
curves denote one - and two-term asymptotic expan- 
sions of the upper boundary for A© 1, k 1. 
(a) A© = 1; (b) Aq = 2/3; (c) A© = 3/5; (d) 
A e = 1/2; (e) A© = 1/4; (f) A© = 1/16. 
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8. Viscous Effects 


The above analysis for nonzero A © ~ 0(e 1/4 ) 
may be extended to include the effects of viscosity. 19 
Adopting the reasonable scalings Pr\ = P ~ 0(1) 
and fi = ~ 0(1) for the liquid-phase Prandtl 

number and the gas-to-liquid viscosity ratio, it was 
previously shown that modifications to the hydro- 
dynamic cellular stability boundary occur, to a first 
approximation, only in the far outer wavenumber 
regime k ~ 0(e _i ), while the leading-order modi- 
fications to the pulsating boundary first occur for 
0(1) wavenumbers. 13,14 Although these results were 
obtained by neglecting the thermal coupling ana- 
lyzed above for the inviscid problem, the same quali- 
tative order-of-magnitude effects remain valid for the 
present case in which A© > 0. Similarly, just as it 
was determined above for the inviscid problem that 
the effects of thermal coupling are felt in the 0(1) 


wavenumber regime for A© = AqC 1 / 4 ~ 0(e 1/M ), a 
similar result is obtained in the viscous case. Thus, 
adopting these scalings, we again seek the neutral 
stability boundary in the 0(1) wavenumber regime, 
anticipating that the primary differences from the 
inviscid analysis presented above will occur on the 
right-hand portion of the stability boundaries ex- 
hibited in Figure 2. In particular, the dispersion 
relation zcj(fc) is sought in the form (65); i.e., as 
iu 0 ~ t~ l / 2 (ivQ + iwxt 1 / 4 4- iu) 2 t 1 ^ 2 F ■ ■ - ). How- 
ever the added algebraic complexity caused by the 
inclusion of the various viscous terms leads us to in- 
troduce corresponding expansions for the coefficients 
b t appearing in Eqs. (41) - (51) directly. 13,14 

The details of the corresponding analysis are 
described elsewhere, 19 but the main result that is 
obtained is that the expression for iu>o is unchanged 
from the inviscid result (70b), and that iu\ is again 
zero. In place of Eq. (71b) for iu 2 , however, we 
obtain the modified expression 

iuo — —2 Pk 2 F -^2 {Ap — P*) 

■ [a; + p* + P ' V 2 k-v 2 A- e (2 a; + pT 3/j ] , 

(91) 

leading to an implicit equation for the viscous pul- 
sating stability boundary A*(fc) given by 

k [(3 + 6)(1 -b) + 8Pfc] 2 = (3 + 6) 2 , 

(92) 

where b - -(2 A*/p* F 1) as before. In the limit 
P — ♦ 0, Eq. (92) collapses to the result (73b). 

Equation (73b) may be transformed to a cu- 
bic equation for fc(6), the single real root of which 
determines the stability boundary. Typical results 
are shown in Figure 5, where, in comparison with 
the corresponding inviscid plot shown in Figure 2, 
it is seen that the main qualitative difference lies 
in the asymptote of the right portion of the bound- 
ary. In particular, that asymptote is no longer given 
simply by b = 1, or equivalently, by A* = -p*, 
but by A* = -p*( 1 F 2PA;) 1/2 , which is the viscous 
pulsating boundary in the limit of zero temperature 
sensitivity. 14 

The pulsating stability boundary (92) may be 
matched to the far outer viscous cellular boundary as 
shown above for the inviscid case, leading to a com- 
posite expansion spanning the outer and far outer 
wavenumber regions. This is exhibited in Figure 6, 
from which we conclude that viscous effects serve to 
widen the stable region for wavenumbers bounded 
away from zero, but do not eliminate the intrinsic 
pulsating instability that exists for sufficiently small 
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Puhntinii Stability Rountlary 
A h ~ 0[t ni ). k - 0(/V- /i' = /, P = 0.0/ 



Figure 5. Pulsating hydrodynamic stability bound- 
aries for k ~ 0(1) and >1© ~ C^e 1 / 4 ). 


Figure 6. Composite hydrodynamic stability bound- 
ary for Ao ~ 0(e 1/f4 ). 


wavenumber disturbances. As in the inviscid case 
(Section 7), the transition from the composite hydro- 
dynamic boundary in Figure 6 to the corresponding 
pulsating and cellular stability boundaries that de- 
scribe the temperature-insensitive limit occurs in the 
regime A© ~ C^e 1 / 2 ) for wavenumbers k ~ 0(e) 1 / 2 . 

9. Summary 

The present work has continued a series of for- 
mal asymptotic treatments of hydrodynamic insta- 
bility in liquid-propellant combustion based on a 
generalized Landau/Levich model that allows for a 
dynamic dependence of the burning rate on local 
pressure and temperature perturbations. The focus 
in the present study was on the effects of thermal 
coupling, represented by nonzero values of an ap- 
propriately defined temperature-sensitivity parame- 
ter Aq, on the nature of a pulsating stability bound- 
ary that exists for negative values of a correspond- 
ing pressure-sensitivity parameter A p . Again ex- 
ploiting the smallness of the gas-to-liquid density 
ratio p, which provides the underlying basis for the 
asymptotic treatment of the stability problem, it was 
shown that for sufficiently large values of A©, the 
pulsating boundary possesses a turning point that 
renders steady, planar burning intrinsically unsta- 
ble for sufficiently small wavenumber perturbations. 


Composite Pulsating! Cellular Stability Boundary 
A p Aq = .5, p* = y = n’ = I , P = .01, e = .005 
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As this parameter decreases, the stable region re- 
emerges such that in the limit that A© vanishes, 
the original separated pulsating and cellular stabil- 
ity branches are recovered. The evolution of the pul- 
sating boundary between these two structures was 
shown to occur on an intermediate wavenumber scale 
relative to the inner, outer and far outer wavenum- 
ber scales that emerge in the asymptotic limit of 
small density ratios in the absence of thermal cou- 
pling. Values of the ratio A©/A p over which this 
evolution occurs were shown to be roughly of the or- 
der of the overall activation energy, suggesting that 
a pulsating, rather than the classical Landau (cel- 
lular) > form of hydrodynamic instability may be the 
more likely manifestation of hydrodynamic instabil- 
ity in at least some types of liquid propellants. 

Acknowledgments 

This work was supported by the United States 
Department of Energy under contract DE-AC04 
-94AL85000 and by the NASA Microgravity Science 
Research Program under contract C-32031-E. 

References 

[1] Landau, L. D., On the Theory of Slow Com- 
bustion, Acta Pkysicochimica URSS 19:77-85 
(1944), Zh . Eksp. i Teor. Fiz. 14:240 (1944). 

[2] Matalon, M., and Matkowsky, B. J., Flames 
as Gasdynamic Discontinuities, J. Fluid Mech. 
124:239-259 (1982). 

[3] Pelce, P., and Clavin, P., Influence of Hydrody- 
namics and Diffusion Upon the Stability Limits 
of Laminar Premixed Flames, J. Fluid Mech. 
124:219-238 (1982). 

[4] Levich, V. G., On the Stability of the Flame 
Front When a Liquid is Burning Slowly, Dokl. 
Akad. Nauk SSSR 109:975-978 (1956). 

[5] Markstein, G. H., Non-Steady Flame Propaga- 
tion , Pergamon Press, Oxford, 1964. 

[6] Sivashinsky, G. I., Diffusional-Thermal Theory 
of Cellular Flames, Combust. Sci. Tech. 15: 
137-146 (1977). 

[7] Margolis, S. B., The Transition to Nonsteady 
Deflagration in Gasless Combustion, Prog. En- 
ergy Combust Sci. 17:135-162 (1991). 

[8] Margolis, S. B., and Williams, F. A., Diffu- 
sional-Thermal Coupling and Intrinsic Instabil- 
ity of Solid Propellant Combustion, Combust 
Sci Tech . 59:27-84 (1988). 


[9] Margolis, S. B., and Williams, F. A., Diffu- 
sional/Thermal Instability of a Solid Propellant 
Flame, SIAM J. AppL Math. 49:1390-1420 
(1989). 

[10] Bechtold, J. K., and Margolis, S. B., Hydrody- 
namic Stability of Solid and Liquid Propellant 
Combustion, Combust Sci Tech. 68:49-84 
(1989). 

[11] Armstrong, R, C., and Margolis, S. B., Hydro- 
dynamic and Reactive/Diffusive Instabilities in 
a Dynamic Model of Liquid Propellant Com- 
bustion, Twenty-Second Symposium (Interna- 
tional) on Combustion , The Combustion Insti- 
tute, Pittsburgh, pp. 1807-1815 (1989). 

[12] Armstrong, R. C., and Margolis, S. B., Hydro- 
dynamic and Reactive/Diffusive Instabilities in 
a Dynamic Model of Liquid Propellant Com- 
bustion — II. Inviscid Fluid Motions, Combust. 
Flame 77:123-138, (1989). 

[13] Margolis, S. B., Hydrodynamic Instability in 
an Extended Landau/Levich Model of Liquid- 
Propellant Combustion at Normal and Reduced 
Gravity, Combust. Flame 113:406-423 (1998). 

[14] Margolis, S. B. } On Pulsating and Cellular 
Forms of Hydrodynamic Instability in Liquid- 
Propellant Combustion, Twenty -Seventh Sym- 
posium (International) on Combustion , pp. 
2375-2386 (1998). 

[15] Vosen, S. R., The Burning Rate of Hydroxy- 
lammonium Nitrate Based Liquid Propellants, 
Twenty-Second Symposium (International) on 
Combustion , The Combustion Institute, Pitts- 
burgh, pp. 1817-1825 (1989). 

[16] Bechtold, J. K., and,Margolis, S. B., Nonlin- 
ear Hydrodynamic Stability and Spinning De- 
flagration of Liquid Propellants, SIAM J . Appl. 
Math . 51:1356-1379 (1991). 

[17] Timnat, Y. M., Advanced Chemical Rocket 
Propulsion , Academic Press, London, 1987. 

[18] Margolis, S. B., Pulsating Hydrodynamic Insta- 
bility and Thermal Coupling in an Extended 
Landau/Levich Model of Liquid-Propellant 
Combustion — I. Inviscid Analysis, to appear 
(1999). 

[19] Margolis, S. B., Pulsating Hydrodynamic Insta- 
bility and Thermal Coupling in an Extended 
Landau/Levich Model of Liquid-Propellant 
Combustion — n. Viscous Analysis, to appear 
(1999). 


19 



